data(faithful)
waiting = as.matrix(faithful[, 2, drop = FALSE])

data(faithful)

waiting  = as.matrix(faithful[, 2, drop = FALSE])

waiting_plot <- ggplot(as.data.frame(waiting), aes(x = waiting)) +
  geom_histogram(aes(y = after_stat(density)), 
                 breaks = seq(min(waiting), max(waiting), length.out = 21),
                 col = "lightcyan4", fill = "lightcyan3") +
  geom_density(col = "indianred2", linewidth = 1
               ) +
  labs(x = "Time between two eruptions (in minutes)", y = "Density") +
  theme_minimal() +
  ggtitle("Density and Histogram of Elapsing Time Between Eruption")

waiting_plot

# prepare the data
X <- waiting
N <- length(X)

# define log.lik helper function
log.lik <- function(X, mu_1, mu_2, var_1, var_2, pi_1, pi_2) {
sum(log(pi_1 * dnorm(X, mu_1, sd = sqrt(var_1)) + pi_2 * dnorm(X, mu_2, sd = sqrt(var_2))))
}

# em function
em_mixture <- function(starting_values, X, tol = .0001, maxits = 100) {
  
  # initialize convergence to false and iteration number to 0
converged <- FALSE 
iter <- 0
N <- length(X)

  # initialize starting values
pi_1 <- starting_values$pi_1 
pi_2 <- starting_values$pi_2 
mu_1 <- starting_values$mu_1 
mu_2 <- starting_values$mu_2 
var_1 <- starting_values$var_1 
var_2 <- starting_values$var_2

  # to save vector of log likelihoods and means
log_likelihoods <- numeric(maxits) 
saved_mu_1 <- numeric(maxits) 
saved_mu_2 <- numeric(maxits) 
ll_changes <- numeric(maxits - 1)
saved_var_1 <- numeric(maxits) 
saved_var_2 <- numeric(maxits) 
saved_pi_1 <- numeric(maxits) 
saved_pi_2 <- numeric(maxits) 
while ((!converged) & (iter < maxits)) {
  
# 1. evaluate the log likelihood at the initial parameters +

  ll <- log.lik(X = X,
                  pi_1 = pi_1,
                  pi_2 = pi_2,
                  mu_1 = mu_1,
                  mu_2 = mu_2,
                  var_1 = var_1,
                  var_2 = var_2)

# 2. E-Step
gamma_1 <- pi_1 * dnorm(X, mu_1, sd = sqrt(var_1)) /
(pi_1 * dnorm(X, mu_1, sd = sqrt(var_1)) + pi_2 * dnorm(X, mu_2, sd = sqrt(var_2)))

gamma_2 <- pi_2 * dnorm(X, mu_2, sd = sqrt(var_2)) /
(pi_1 * dnorm(X, mu_1, sd = sqrt(var_1)) + pi_2 * dnorm(X, mu_2, sd = sqrt(var_2)))

# 3. M-Step
# 
pi_1 <- sum(gamma_1)/N 
pi_2 <- sum(gamma_2)/N
mu_1 <- sum(X * gamma_1) / sum(gamma_1)
mu_2 <- sum(X * gamma_2) / sum(gamma_2)
var_1 <- sum((X - mu_1)^2 * gamma_1) / sum(gamma_1) 
var_2 <- sum((X - mu_2)^2 * gamma_2) / sum(gamma_2)

    # 4. evaluate the log likelihood at the new parameter values
ll.new <- log.lik(X = X,
pi_1 = pi_1,
                      pi_2 = pi_2,
                      mu_1 = mu_1,
                      mu_2 = mu_2,
                      var_1 = var_1,
                      var_2 = var_2)

    # store new parameters
log_likelihoods[iter + 1] <- ll.new 
saved_mu_1[iter + 1] <- mu_1 
saved_mu_2[iter + 1] <- mu_2
saved_var_1[iter + 1] <- var_1
saved_var_2[iter + 1] <- var_2
saved_pi_1[iter + 1] <- pi_1
saved_pi_2[iter + 1] <- pi_2

if (iter > 0) {
ll_changes[iter] <- abs(ll - ll.new) }

    # 5. check convergence
if(abs(ll - ll.new) < tol) { converged <- TRUE
}

    # next iteration
iter <- iter + 1

    # message to keep track of progress
cat(paste0("Running iteration ", iter, ". Log likelihood changed by ", round(abs(ll - ll.new), 4), "\n"))
}

  # save the parameter values
params <- data.frame(pi_1 = pi_1, 
                     pi_2 = pi_2,
                     mu_1 = mu_1,
                     mu_2 = mu_2,
                     var_1 = var_1, 
                     var_2 = var_2) %>%
  pivot_longer(
    cols = everything(),
    names_to = c(".value", "group"), names_pattern = "(.*)_(.)"
    )

final_gammas_df <- data.frame( obs = 1:length(gamma_1),
waiting_1 = round(gamma_1, 5),
waiting_2 = round(gamma_2, 5) )

return(list(
params = params,
log_likelihoods = log_likelihoods[1:iter], 
ll_changes = ll_changes[1:(iter - 1)], 
saved_mu_1 = saved_mu_1[1:iter], 
saved_mu_2 = saved_mu_2[1:iter], 
final_gammas = (head(final_gammas_df)),
saved_var_1 = saved_var_1[1:iter],
saved_var_2 = saved_var_2[1:iter],
saved_pi_1 = saved_pi_1[1:iter],
saved_pi_2 = saved_pi_2[1:iter]
))
}
starting_values_1 <- list(pi_1 = .5, pi_2 = .5, mu_1 = 55, mu_2 = 80, var_1 = 1, var_2 = 1)
starting_values_2 <- list(pi_1 = .5, pi_2 = .5, mu_1 = 30, mu_2 = 60, var_1 = 1, var_2 = 1)

em1 <- em_mixture(starting_values = starting_values_1, X = waiting)
## Running iteration 1. Log likelihood changed by 3842.1989
## Running iteration 2. Log likelihood changed by 0.2425
## Running iteration 3. Log likelihood changed by 0.0265
## Running iteration 4. Log likelihood changed by 0.0101
## Running iteration 5. Log likelihood changed by 0.0043
## Running iteration 6. Log likelihood changed by 0.0019
## Running iteration 7. Log likelihood changed by 0.0008
## Running iteration 8. Log likelihood changed by 0.0003
## Running iteration 9. Log likelihood changed by 0.0001
## Running iteration 10. Log likelihood changed by 0.0001
em1
## $params
## # A tibble: 2 × 4
##   group    pi    mu   var
##   <chr> <dbl> <dbl> <dbl>
## 1 1     0.361  54.6  34.5
## 2 2     0.639  80.1  34.4
## 
## $log_likelihoods
##  [1] -1034.288 -1034.046 -1034.019 -1034.009 -1034.005 -1034.003 -1034.002
##  [8] -1034.002 -1034.002 -1034.002
## 
## $ll_changes
## [1] 0.24249928394 0.02654319958 0.01006483672 0.00430799398 0.00185608116
## [6] 0.00080090007 0.00034596748 0.00014956708 0.00006469558
## 
## $saved_mu_1
##  [1] 54.75000 54.73771 54.69984 54.67111 54.65179 54.63909 54.63077 54.62531
##  [9] 54.62173 54.61938
## 
## $saved_mu_2
##  [1] 80.28488 80.18128 80.14547 80.12625 80.11417 80.10628 80.10109 80.09767
##  [9] 80.09541 80.09393
## 
## $final_gammas
##   obs waiting waiting.1
## 1   1 0.00011   0.99989
## 2   2 0.99991   0.00009
## 3   3 0.00420   0.99580
## 4   4 0.96774   0.03226
## 5   5 0.00000   1.00000
## 6   6 0.99981   0.00019
## 
## $saved_var_1
##  [1] 34.40750 35.54276 35.31765 35.04204 34.84594 34.71642 34.63183 34.57657
##  [9] 34.54040 34.51668
## 
## $saved_var_2
##  [1] 31.48280 33.29762 33.79021 34.02135 34.16164 34.25301 34.31332 34.35318
##  [9] 34.37949 34.39684
## 
## $saved_pi_1
##  [1] 0.3676471 0.3648946 0.3634577 0.3625668 0.3619894 0.3616113 0.3613630
##  [8] 0.3611997 0.3610924 0.3610218
## 
## $saved_pi_2
##  [1] 0.6323529 0.6351054 0.6365423 0.6374332 0.6380106 0.6383887 0.6386370
##  [8] 0.6388003 0.6389076 0.6389782
em2 <- em_mixture(starting_values = starting_values_2, X = waiting)
## Running iteration 1. Log likelihood changed by 40474.0626
## Running iteration 2. Log likelihood changed by 1.4258
## Running iteration 3. Log likelihood changed by 1.5858
## Running iteration 4. Log likelihood changed by 1.3444
## Running iteration 5. Log likelihood changed by 0.6242
## Running iteration 6. Log likelihood changed by 0.287
## Running iteration 7. Log likelihood changed by 0.1564
## Running iteration 8. Log likelihood changed by 0.105
## Running iteration 9. Log likelihood changed by 0.084
## Running iteration 10. Log likelihood changed by 0.0753
## Running iteration 11. Log likelihood changed by 0.0736
## Running iteration 12. Log likelihood changed by 0.0785
## Running iteration 13. Log likelihood changed by 0.0917
## Running iteration 14. Log likelihood changed by 0.1174
## Running iteration 15. Log likelihood changed by 0.1616
## Running iteration 16. Log likelihood changed by 0.2321
## Running iteration 17. Log likelihood changed by 0.3354
## Running iteration 18. Log likelihood changed by 0.4769
## Running iteration 19. Log likelihood changed by 0.6618
## Running iteration 20. Log likelihood changed by 0.8938
## Running iteration 21. Log likelihood changed by 1.174
## Running iteration 22. Log likelihood changed by 1.5203
## Running iteration 23. Log likelihood changed by 1.9993
## Running iteration 24. Log likelihood changed by 2.6826
## Running iteration 25. Log likelihood changed by 3.4622
## Running iteration 26. Log likelihood changed by 4.0058
## Running iteration 27. Log likelihood changed by 4.1515
## Running iteration 28. Log likelihood changed by 4.0788
## Running iteration 29. Log likelihood changed by 3.9936
## Running iteration 30. Log likelihood changed by 3.9807
## Running iteration 31. Log likelihood changed by 4.0246
## Running iteration 32. Log likelihood changed by 3.9827
## Running iteration 33. Log likelihood changed by 3.6752
## Running iteration 34. Log likelihood changed by 3.0912
## Running iteration 35. Log likelihood changed by 2.399
## Running iteration 36. Log likelihood changed by 1.7363
## Running iteration 37. Log likelihood changed by 1.1553
## Running iteration 38. Log likelihood changed by 0.6942
## Running iteration 39. Log likelihood changed by 0.3767
## Running iteration 40. Log likelihood changed by 0.188
## Running iteration 41. Log likelihood changed by 0.0886
## Running iteration 42. Log likelihood changed by 0.0403
## Running iteration 43. Log likelihood changed by 0.018
## Running iteration 44. Log likelihood changed by 0.0079
## Running iteration 45. Log likelihood changed by 0.0035
## Running iteration 46. Log likelihood changed by 0.0015
## Running iteration 47. Log likelihood changed by 0.0007
## Running iteration 48. Log likelihood changed by 0.0003
## Running iteration 49. Log likelihood changed by 0.0001
## Running iteration 50. Log likelihood changed by 0.0001
em2
## $params
## # A tibble: 2 × 4
##   group    pi    mu   var
##   <chr> <dbl> <dbl> <dbl>
## 1 1     0.361  54.6  34.4
## 2 2     0.639  80.1  34.5
## 
## $log_likelihoods
##  [1] -1095.345 -1093.919 -1092.334 -1090.989 -1090.365 -1090.078 -1089.922
##  [8] -1089.817 -1089.733 -1089.657 -1089.584 -1089.505 -1089.414 -1089.296
## [15] -1089.135 -1088.902 -1088.567 -1088.090 -1087.428 -1086.535 -1085.360
## [22] -1083.840 -1081.841 -1079.158 -1075.696 -1071.690 -1067.539 -1063.460
## [29] -1059.466 -1055.486 -1051.461 -1047.478 -1043.803 -1040.712 -1038.313
## [36] -1036.577 -1035.421 -1034.727 -1034.351 -1034.163 -1034.074 -1034.034
## [43] -1034.016 -1034.008 -1034.004 -1034.003 -1034.002 -1034.002 -1034.002
## [50] -1034.002
## 
## $ll_changes
##  [1] 1.42580723273 1.58578318422 1.34440088607 0.62422788506 0.28697213883
##  [6] 0.15638497314 0.10501440105 0.08401235063 0.07528822039 0.07361061792
## [11] 0.07847914200 0.09174982788 0.11739559141 0.16164337823 0.23208819827
## [16] 0.33544312382 0.47688142788 0.66181112993 0.89381152691 1.17400170099
## [21] 1.52029439942 1.99928550981 2.68262673657 3.46222489501 4.00583768251
## [26] 4.15148012708 4.07878500674 3.99355798650 3.98070687911 4.02461514856
## [31] 3.98265459756 3.67520501228 3.09121240806 2.39895017144 1.73627559810
## [36] 1.15530292225 0.69416182073 0.37665394614 0.18803046515 0.08860539929
## [41] 0.04029663788 0.01795439269 0.00790726233 0.00345934260 0.00150753994
## [46] 0.00065542944 0.00028454811 0.00012342167 0.00005350274
## 
## $saved_mu_1
##  [1] 44.20000 44.99270 45.57248 45.84219 45.95954 46.03948 46.10221 46.15491
##  [9] 46.20275 46.24875 46.29481 46.34275 46.39498 46.45487 46.52708 46.61744
## [17] 46.73233 46.87780 47.05950 47.28271 47.55098 47.86586 48.23165 48.65997
## [25] 49.15371 49.67466 50.15822 50.57768 50.95876 51.34081 51.74386 52.16116
## [33] 52.57092 52.95541 53.30660 53.61766 53.88038 54.09020 54.24919 54.36463
## [41] 54.44583 54.50168 54.53952 54.56491 54.58183 54.59305 54.60048 54.60538
## [49] 54.60861 54.61075
## 
## $saved_mu_2
##  [1] 71.14471 71.10225 71.17822 71.32847 71.47849 71.58064 71.64511 71.68783
##  [9] 71.71970 71.74654 71.77130 71.79569 71.82121 71.84957 71.88314 71.92512
## [17] 71.97950 72.05078 72.14393 72.26463 72.41922 72.61440 72.85895 73.16941
## [25] 73.57151 74.08209 74.68207 75.32268 75.95841 76.56663 77.14541 77.69499
## [33] 78.20031 78.63903 79.00268 79.29692 79.52894 79.70443 79.83091 79.91840
## [41] 79.97727 80.01628 80.04195 80.05881 80.06987 80.07714 80.08191 80.08504
## [49] 80.08711 80.08846
## 
## $final_gammas
##   obs waiting waiting.1
## 1   1 0.00010   0.99990
## 2   2 0.99991   0.00009
## 3   3 0.00408   0.99592
## 4   4 0.96705   0.03295
## 5   5 0.00000   1.00000
## 6   6 0.99981   0.00019
## 
## $saved_var_1
##  [1]  0.9600000  1.0268115  0.7063628  0.5077563  0.4906262  0.5289275
##  [7]  0.5806714  0.6373783  0.6972916  0.7605596  0.8287625  0.9049071
## [13]  0.9936607  1.1017980  1.2386355  1.4159806  1.6473870  1.9477409
## [19]  2.3335815  2.8214029  3.4255237  4.1654256  5.0841950  6.2334865
## [25]  7.5842371  8.9823096 10.3210570 11.6557531 13.0816774 14.6591372
## [31] 16.4057082 18.2822749 20.2397722 22.2756492 24.3720843 26.4245451
## [37] 28.2933362 29.8826164 31.1597204 32.1386510 32.8598983 33.3745051
## [43] 33.7327431 33.9776927 34.1430870 34.2538084 34.3275022 34.3763624
## [49] 34.4086749 34.4300077
## 
## $saved_var_2
##  [1] 179.17015 180.23681 178.98106 176.31071 173.58818 171.73370 170.57450
##  [8] 169.81798 169.26369 168.80405 168.38538 167.97704 167.55413 167.08822
## [15] 166.54147 165.86308 164.99006 163.85261 162.37669 160.48010 158.07309
## [22] 155.06235 151.32685 146.62998 140.55652 132.69332 123.04488 112.20023
## [29] 101.03600  90.25977  80.14728  70.74596  62.28422  55.14939  49.48332
## [36]  45.11197  41.80520  39.38798  37.69742  36.56008  35.81374  35.32951
## [43]  35.01611  34.81290  34.68073  34.59448  34.53807  34.50110  34.47684
## [50]  34.46091
## 
## $saved_pi_1
##  [1] 0.009191176 0.007858950 0.010980381 0.016927210 0.022784483 0.026763999
##  [7] 0.029285957 0.030970839 0.032238956 0.033315887 0.034315485 0.035305789
## [13] 0.036346204 0.037508260 0.038889551 0.040622609 0.042873693 0.045831816
## [19] 0.049707049 0.054742493 0.061209144 0.069391696 0.079663450 0.092713280
## [25] 0.109528875 0.130494476 0.154340009 0.178849101 0.202456914 0.224752634
## [31] 0.245982961 0.266232466 0.284956194 0.301436265 0.315441921 0.327106898
## [37] 0.336544468 0.343846669 0.349227855 0.353033569 0.355648036 0.357411789
## [43] 0.358589106 0.359370286 0.359886876 0.360227834 0.360452618 0.360600711
## [49] 0.360698236 0.360762444
## 
## $saved_pi_2
##  [1] 0.9908088 0.9921411 0.9890196 0.9830728 0.9772155 0.9732360 0.9707140
##  [8] 0.9690292 0.9677610 0.9666841 0.9656845 0.9646942 0.9636538 0.9624917
## [15] 0.9611104 0.9593774 0.9571263 0.9541682 0.9502930 0.9452575 0.9387909
## [22] 0.9306083 0.9203366 0.9072867 0.8904711 0.8695055 0.8456600 0.8211509
## [29] 0.7975431 0.7752474 0.7540170 0.7337675 0.7150438 0.6985637 0.6845581
## [36] 0.6728931 0.6634555 0.6561533 0.6507721 0.6469664 0.6443520 0.6425882
## [43] 0.6414109 0.6406297 0.6401131 0.6397722 0.6395474 0.6393993 0.6393018
## [50] 0.6392376


Indeed, for EM2

# for em2
em_params_2 <- em2$params
em_log_likelihoods_2 <- em2$log_likelihoods
em_saved_mu_1_2 <- em2$saved_mu_1
em_saved_mu_2_2 <- em2$saved_mu_2
em_final_gammas_2 <- em2$final_gammas
ll_changes_2 <- em2$ll_changes

# plotting log likelihoods
em_ll_df_2 <-  tibble(
  obs = 1:length(em_log_likelihoods_2),
  em_log_likelihoods = (em_log_likelihoods_2)
 )

em_ll_plot_2 <- ggplot(em_ll_df_2, aes(x = obs, y = em_log_likelihoods)) +
  geom_line(col = "indianred2", linewidth = 1) +  
  theme_minimal()

em_ll_plot_2

# plotting change in log likelihood at each iteration
ll_change_df_2 <- data.frame(
    iteration = 1:length(ll_changes_2),
    ll_change = ll_changes_2
)

ggplot(ll_change_df_2, aes(x = iteration, y = ll_change)) +
    geom_line() +
    theme_minimal()


While for EM1:

# for em1
em_params_1 <- em1$params
em_log_likelihoods_1 <- em1$log_likelihoods
em_saved_mu_1_1 <- em1$saved_mu_1
em_saved_mu_2_1 <- em1$saved_mu_2
em_final_gammas_1 <- em1$final_gammas
ll_changes_1 <- em1$ll_changes

# plotting log likelihoods
em_ll_df_1 <-  tibble(
  iteration = 1:length(em_log_likelihoods_1),
  em_log_likelihoods = (em_log_likelihoods_1)
 )

em_ll_plot_1 <- ggplot(em_ll_df_1, aes(x = iteration, y = em_log_likelihoods)) +
  geom_line(col = "indianred2", linewidth = 1) +  
  theme_minimal()

em_ll_plot_1

# plotting change in log likelihood at each iteration
ll_change_df_1 <- data.frame(
    iteration = 1:length(ll_changes_1),
    ll_change = ll_changes_1
)

ggplot(ll_change_df_1, aes(x = iteration, y = ll_change)) +
    geom_line() +
    theme_minimal()


Asya’s Fun

It plots the log likelihood as a function of mu_1 at a particular iteration (evaluated at the current maximizing levels of the other parameters)

as a function of mu_1

log.lik2 <- function(data, mu_1, mu_2, var_1, var_2, pi_1, pi_2) {
  sum(log(pi_1 * dnorm(data, mu_1, sd = sqrt(var_1)) + pi_2 * dnorm(data, mu_2, sd = sqrt(var_2))))
}

plot.lik.tot.mu_1 <- function(iter, em) {
  x_values <- seq(0, 95, length.out = 1000)
  ll <- sapply(x_values, function(x) log.lik2(data = waiting, 
                                              mu_1 = x,
                                              mu_2 = em$saved_mu_2[iter],
                                              var_1 = em$saved_var_1[iter], 
                                              var_2 = em$saved_var_2[iter],
                                              pi_1 = em$saved_pi_1[iter],
                                              pi_2 = em$saved_pi_2[iter]))
  plot(x_values, ll, type = "l")
  points(x_values[which.max(ll)], max(ll), col = "red", pch = 19)
}

plot.lik.part.mu_1 <- function(iter, em) {
  ll <- sapply(X = seq(40, 60, length.out = 1000),
               FUN = function(x) log.lik2(data = waiting,
                                         mu_1 = x,
                                         mu_2 = em$saved_mu_2[iter],
                                         var_1 = em$saved_var_1[iter],
                                         var_2 = em$saved_var_2[iter],
                                         pi_1 = em$saved_pi_1[iter],
                                         pi_2 = em$saved_pi_2[iter]))
  
  plot(seq(40, 60, length.out = 1000), ll, type = "l")
  points(seq(40, 60, length.out = 1000)[which.max(ll)], max(ll), col = "red", pch = 19)
}

plot.lik.tot.mu_1(em = em2, iter = 9)

plot.lik.part.mu_1(em = em2, iter = 9)

as a function of mu_2

plot.lik.tot_mu2 <- function(iter, em) {
  x_values <- seq(0, 95, length.out = 1000)
  ll <- sapply(x_values, function(x) log.lik2(data = waiting, 
                                              mu_1 = em$saved_mu_1[iter],
                                              mu_2 = x, 
                                              var_1 = em$saved_var_1[iter], 
                                              var_2 = em$saved_var_2[iter],
                                              pi_1 = em$saved_pi_1[iter],
                                              pi_2 = em$saved_pi_2[iter]))
  plot(x_values, ll, type = "l")
  points(x_values[which.max(ll)], max(ll), col = "red", pch = 19)
}

plot.lik.part_mu2 <- function(iter, em) {
  ll <- sapply(X = seq(40, 60, length.out = 1000),
               FUN = function(x) log.lik2(data = waiting,
                                         mu_1 = em$saved_mu_1[iter],
                                         mu_2 = x,
                                         var_1 = em$saved_var_1[iter],
                                         var_2 = em$saved_var_2[iter],
                                         pi_1 = em$saved_pi_1[iter],
                                         pi_2 = em$saved_pi_2[iter]))
  
  plot(seq(40, 60, length.out = 1000), ll, type = "l")
  points(seq(40, 60, length.out = 1000)[which.max(ll)], max(ll), col = "red", pch = 19)
}

plot.lik.tot_mu2(em = em2, iter = 14)

plot.lik.part_mu2(em = em2, iter = 14)

Given fixed parameters

plot.lik.tot <- function(iter, em) {
  x_values <- seq(0, 95, length.out = 1000)
  ll <- sapply(x_values, function(x) log.lik2(data = x, 
                                              mu_1 = em$saved_mu_1[iter],
                                              mu_2 = em$saved_mu_2[iter],
                                              var_1 = em$saved_var_1[iter], 
                                              var_2 = em$saved_var_2[iter],
                                              pi_1 = em$saved_pi_1[iter],
                                              pi_2 = em$saved_pi_2[iter]))
  plot(x_values, ll, type = "l")
  points(x_values[which.max(ll)], max(ll), col = "red", pch = 19)
}

plot.lik.part <- function(iter, em) {
  ll <- sapply(X = seq(40, 60, length.out = 1000),
               FUN = function(x) log.lik2(data = x,
                                         mu_1 = em$saved_mu_1[iter],
                                         mu_2 = em$saved_mu_2[iter],
                                         var_1 = em$saved_var_1[iter],
                                         var_2 = em$saved_var_2[iter],
                                         pi_1 = em$saved_pi_1[iter],
                                         pi_2 = em$saved_pi_2[iter]))
  
  plot(seq(40, 60, length.out = 1000), ll, type = "l")
  points(seq(40, 60, length.out = 1000)[which.max(ll)], max(ll), col = "red", pch = 19)
}

plot.lik.tot(em = em2, iter = 14)

plot.lik.part(em = em2, iter = 14)

total animation EM2 fixed

# complete plot
plot.lik.tot.animation <- function(em, start_iter, end_iter) {
  x_values <- seq(0, 90, length.out = 1000)
  data_list <- lapply(start_iter:end_iter, function(iter) {
    ll <- sapply(x_values, FUN = function(x) log.lik2(data = x,
                                                      mu_1 = em$saved_mu_1[iter],
                                                      mu_2 = em$saved_mu_2[iter],
                                                      var_1 = em$saved_var_1[iter],
                                                      var_2 = em$saved_var_2[iter],
                                                      pi_1 = em$saved_pi_1[iter],
                                                      pi_2 = em$saved_pi_2[iter]))
    data.frame(iter = iter, x = x_values, ll = ll)
  })

  data <- do.call(rbind, data_list)

  p <- ggplot(data, aes(x = x, y = ll, group = iter)) +
       geom_line() +
       transition_states(states = iter, transition_length = 2, state_length = 1) +
       labs(title = "Iteration: {closest_state}") +
       theme_minimal()

  animate(p, nframes = 100, duration = 5, width = 800, height = 600)
}

# Use the function like this
plot.lik.tot.animation(em = em2, start_iter = 0, end_iter = 50)

total animation EM1 fixed

plot.lik.tot.animation(em = em1, start_iter = 0, end_iter = 10)

total animation EM2 varying mu_1

# complete plot
plot.lik.tot.mu_1_animation <- function(em, start_iter, end_iter) {
  x_values <- seq(0, 95, length.out = 1000)
  data_list <- lapply(start_iter:end_iter, function(iter) {
    ll <- sapply(x_values, function(x) log.lik2(data = waiting, 
                                                mu_1 = x,
                                                mu_2 = em$saved_mu_2[iter],
                                                var_1 = em$saved_var_1[iter], 
                                                var_2 = em$saved_var_2[iter],
                                                pi_1 = em$saved_pi_1[iter],
                                                pi_2 = em$saved_pi_2[iter]))
    data.frame(iter = iter, x = x_values, ll = ll)
  })

  data <- do.call(rbind, data_list)
  p <- ggplot(data, aes(x = x, y = ll, group = iter)) +
       geom_line() +
       transition_states(states = iter, transition_length = 2, state_length = 1) +
       labs(title = "Iteration: {closest_state}") +
       theme_minimal()

  animate(p, duration = 5, nframes = 100, width = 800, height = 600)
}

# Use the function like this
plot.lik.tot.mu_1_animation(em = em2, start_iter = 1, end_iter = 30)